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ABSTRACT 

Aims. Theory of random processes provides an attractive mathematical tool to describe the fluctuating signal from accreting sources, such as 
active galactic nuclei and Galactic black holes observed in X-rays. These objects exhibit featureless variability on different timescales, probably 
originating from an accretion disc. 

Methods. We study the basic features of the power spectra in terms of a general framework, which permits semi-analytical determination of 
the power spectral density (PSD) of the resulting light curve. We consider the expected signal generated by an ensemble of spots randomly 
created on the accretion disc surface. Spot generation is governed by Poisson or by Hawkes processes. The latter one represents an avalanche 
mechanism and seems to be suggested by the observed form of the power spectrum. We include general relativity effects shaping the signal on 
its propagation to a distant observer 

Results. We analyse the PSD of a spotted disc light curve and show the accuracy of our semi-analytical approach by comparing the obtained 
PSD with the results of Monte Carlo simulations. The asymptotic slopes of PSD are at low frequencies and they drop to -2 at high frequencies, 
usually with a single frequency break. More complex two-peak solutions also occur. The amplitude of the peaks and their frequency difference 
depend on the inherent timescales of the model, i.e., the intrinsic lifetime of the spots and the typical duration of avalanches. 
Conclusions. At intermediate frequencies, the intrinsic PSD is influenced by the individual light curve profile as well as by the type of the 
underlying process. However, even in cases when two Lorentzians seem to dominate the PSD, it does not necessarily imply that two single 
oscillation mechanisms operate simultaneously. Instead, it may well be the manifestation of the avalanche mechanism. The main advantage of 
our approach is an insight in the model functioning and the fast evaluation of the PSD. 
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1. Introduction 

It is widely accepted that massive black holes with accretion 
discs reside in cores of active galactic nuclei, where most ac- 



tivity originates and X-rays are produced (e.g.. lBlandford et al 



11990.) . The observed light curves, / = f(t), show irregu- 
lar, featureless fluctuations with a very complex behaviour , 
practically at every studied frequency f Gaskell et al.l 12006). 
Variabili ty has been traditionall y analysed by the Fourier 
method dFeigelson & Babulll992l) . Remarkably, a number of 
similarities appear between the properties of massive black 
holes in galactic nuclei and those in X-ray binaries, suggesting 
that some kind of universal r escaling operates according t o cen- 
tral masses of these systems dMirabel & Rodrig uez 1998) . This 



the 'rms' characteristic), function of a single variable (for ex- 
ample, the power-spectral density - PSD) or of many vari- 
ables (e.g., poly-spectra, rms-flux relation, etc). The appropri- 
ate choice depends on the type of information we seek and the 
quality of data available. The PSD is a traditional and widely 
utilised method to examine variable signals, and the AGN light 
curves are no exception. A typical signal can be represented 
by a broad band PSD with the tendency towards flattening 



concerns also the X-ray po wer spectra (e.g., iMarkowitz et al 
2003l:lMcHardv et alJl2006l) . 



Light curves can be characterised by an appropriate es- 
timator of the source variability which, in the mathematical 
sense, is a functional: / — > S[f]. We accept the idea that 
the signal resulting from a spotted accretion disc is intrin- 
sically stochastic, likely originating from turbulence. Hence, 
S[f] is a random value. It can be a number (for example, 



at low frequency f Lawrence et al. 198 T . McHardv & Czerny 
1987; Lawrence & Pap adakisrfr993l: iMushotzkv et alj Il993[ 
lUttlev et alJl2002h . 



There is an ongoing debate about the overall shape of the 
PSD and the occurrence of the break frequency or, possibly, 
two break frequenci e s at which the slope of th e PSD can change 
dNowak et alJll999l: IMarkowitz et al.ll2003h . In the case of a 



widely -studied Seyfert galaxy, MCG-6-30-15. lMcHardy et al 



(l2005h have closely examined the slope of PSD, namely its 
bending, with RXTE and XMM-Newton data. It is worth notic- 
ing that the accurate fits to the X-ray sources seem to exhibit a 
muIti-Lorentzian structure rather than a simple power law. The 
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same is tr ue for the best studied example, the Seyfert 1 galaxy 

Akn 564 dMcHardv et al.ll2007h. 

It was pr oposed (^ Abramowicz et alj Il99ll : IZhang & Bad 

I991I: IWiita e t al. 1992|) that hot spots contribute to the vari- 



abihty of the AGN variability, or that they could even be 
the dominant process shaping the variability pattern. These 
spots should occur on the disc surface following its inter- 
mittent irradiation by localised coronal flares (iGaleev et alJ 
19791; iMerloni & FabianlbOOll; ICzernv et alj|2004 . Here, the 
"spots" represent a somewhat generic designation for non- 
axisymmetric features evolving on the disc surface in connec- 
tion with flares. They share the bulk orbital motion with the 
underlying disc. The observed signal is thus modulated by rela- 
tivistic effects as photons propagate towards a distant observer. 
Various schemes have been discussed in which the fluctu- 
ations of the disc emissivity at different points of space and 
time are mutually inte rconnected in some way. In particular , 
the avalanche model (Pout anen & Fabianlll999l: IZyckil l2002t 
IZyck i & Niedzwiecki 2005) seems to be physically substan- 
tiated within the framework of magnetically-triggered flares 
and spots. It is also a promising model capable to repro- 
duce, for example, the broken power-law PSD profiles. Notice, 
h owever, that other promising ideas were proposed (e.g.. 



Mineshige et al.lll994c lLyubarskiilll997h . provoking the ques 



tion of whether a common mathematical basis could reflect the 
entire range of models and provide us with general constraints, 
independent of (largely unknown) model details. 

We add to th is model by applyin g the method of random 
point processes (ICox & Miller! 1 1965h . Interestingly enough, a 
rather formal approach can provide useful analytical formu- 
lae defining the basic form of the expected power spectrum. 
Apart from this practical aspect, we suggest that the concept 
discussed here offers much better insight into various influ- 
ences that shape the expected form the power spectrum. These 
are very attractive features especially with respect to avalanche 
models, which may have different flavours, typically with a vast 
parameter space, thus proving very demanding to examine in a 
systematic manner 

Even more important is that the adopted formahsm pro- 
vides a very general tool and allows for a broader perspec- 
tive o n different mechanisms of variability dPechacek & Karas 
20071). We develop the idea in a systematic way and 



give the explicit coiTespondence between our approach and 
some of the above-men t ioned and widely-k nown scenarios 
dAbramowicz et al.lll99ll : IPoutanen & Fabianlll999;) . This de- 
scription provides semi-analytical solutions, convenient to 
search through a broad parameter space. Our results can help to 
identify how the intrinsic properties of individual flares and the 
relativistic effects influence the overall PSD. In particular, we 
can identify those situations in which a doubly-broken power 
law occurs. 

We consider stochastic processes (e.g 



FeUer 1971 



Gardinenl994 ') in the framework appropriate for modelling the 



accretion disc variability. In particular, in Sect. |2] we consider 
a simple model of a spotted accretion disc constrained by the 
following three assumptions about the creation and evolution 
of spots; (i) each spot is described by its time and place of birth 
{tj, rj and (pj) in the plain of the disc; (ii) every new occurrence 



starts instantaneously; afterwards, the emissivity decays grad- 
ually to zero (the total radiated energy is of course finite); and 
(iii) the intrinsic emissivity is fully determined by a finite set 
of parameters which form a vector, ^ , defining the light curve 
profile. Later on, we will consider the modulation of the in- 
trinsic emission by the orbital motion and relativistic lensing. 
The disc itself has a passive role in our considerations; we will 
treat it as a geometrically thin, optically thick layer lying in the 
equatorial plane. 

Because of the apparently random nature of the variabil- 
ity, we adopt a stochastic model in which the creation of spots 
is governed by a random process. The assumption that spots 
are mutually statistically independent seems to be a reasonable 
(first) approximation, however, we find that we do need to in- 
troduce some kind of relationship between them. This connec- 
tion is discussed in Sect. [3] The statistical dependence among 
spots can be introduced in several ways. In Sect. H] we ex- 
plore in detail the specific models of interrelated spots using 
the formalism of Hawkes-type processes. Conclusions are sum- 
marised in Sect. |5] Finally, in the Appendix we provide some 
mathematical prerequisites, which the reader may find useful 
to understand the general background of the paper, and we also 
summarise the mathematical notation. 

2. Models of stochastic variability 

2.1. Orbiting spots and relativistic effects 

We will apply our investigations to models where the signal is 
produced by point-like orbiting spots (circular Keplerian mo- 
tion along the azimuthal (/> direction). The intrinsic emission, 
produced in the local co-orbiting frame of the spot, is influ- 
enced by the Doppler effect and gravitational lensing, which 
cannot be ignored at typical distances of several units or tens 
of gravitational radii. Photons emitted at different moments and 
positions experience different light-travel time on the way to- 
wards the observer, so the observed timing properties should 
reflect this specific modulation. We adopt the Schwarzschild 
metric for the gr avitational field and employ the method of 
transfer function (IPechacek et al.l 120051 120061) to describe the 
light amplification (or dilution); 60 is inclination angle of the 
observer {9o - 90 deg corresponds the edge-on view of the 
disc plane). The periodical modulation of the observed signal 
is included in the transfer function F{t, r, 0o) = F{(p{t), r, 0o). 
An implicit relation holds for the phase. 



<^(0 = Qf - 6t(<l){t), r, 0o), 



(1) 



where the effect of time delays 5f(0, r, 0o) is taken into ac- 
count. The modulation by F is superimposed onto that due 
to intrinsic timescales present in the signal from spots. Then, 
for the final flux f{t) measured by a distant observer, we write 
fit) ^ IF{t,r,B,). 

In the case of an infinitesimal surface element with in- 
trinsically constant and isotropic emissivity /, orbiting with 
Keplerian orbital frequency 0.{r), the flux measured by a dis- 
tant observer varies just as F changes along the orbit. 

We remind the reader that the mass of central black holes 
in galactic nuclei is in the range of ^ M » IO^-IO'Mq. 
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Mass of the accretion disc is at least three orders of mag- 
nitude smaller than the black hole mass, so we neglect it in 
our calculations (the accretion disc self-gravity may be im- 
portant for its intrinsic struct ure, but the direct gravitational 
effect on light is quite small; iKaras et al.lll995h . Hence, the 
gravitational fie ld can indeed be de scribed by a vacuum black- 
hole spacetime (IMisner et al.lll973l) . We use geometrical units 
with c = G = 1. Transformation to physical scales can be 
achieved when the mass of the central black hole is specified 
because Keplerian frequency scales inversely with the gravita- 
tional radius. The gravitational radius of a massive black hole 
is Rg = c'^GM ^ 0.48 x lO^^Mgpc, and the corresponding 
characteristic time-scale is fg = c^^GM » 0.49 x lO^Ms sec, 
where the mass Mg = M/{IQ^Mq). All lengths and times can 
be made dimensionless by expressing them in units of M, so 
they can be easily scaled to different masses. For example, for 
the Keplerian orbital period we obtain To^b ~ 3.1 x lO"* P^^M^, 
where the radius is expressed in units of R^ and Toib is given in 
seconds. 

Let us note that the intrinsic timescales of the spot evolution 
and of avalanches (both timescales will be discussed below) 
need not to be directly connected with the Keplerian orbital pe- 
riod. This internal freedom of the model can help to bring the 
predicted frequency of the breaks of the PSD profile in agree- 
ment with the data. 



2.2. Model driven by a general point process 

Now we will describe the process of the creation of spots from 
the statistics point of view. Let us consider a signal of the form 



/(r) = 2/(f-<J„^pF(f-<Jp,-,o). 



(2) 



The underlying process consists of a sequence of events which, 
in general, can be either mutually independent, or there can 
be some statistical dependence among them. Naturally, the lat- 
ter case will be more complicated and interesting. In Eq. (|2|i, 
I(t, ^) - g(t, ^) 6{t) is the light curve profile of a single event; 
6j - tj + tdj and 6pj - 6j + t^j are the time offsets (determined by 
the moment of the ignition of the spot) and the initial time de- 
lay (which is an arbitrary but fixed value); 0{t) is the Heaviside 
function; and g(t, f ) is a non-negative function of k + I vari- 
ables, t and ^ - (^\ . . . ,^'^). Hereafter, we omit the explicit 
dependency on 6o. 

Quantities ^ , tj, rj, tpj and tdj are random values. The 
vector ^ determines the duration and shape of each event (tj 
is time of ignition; parameter tpj determines the initial phase 
of the periodical modulation of the 7-th event; and t^j is the 
corresponding initial time-offset). These assumptions bring the 
formulation of the problem close to the framework studied by 
Bremaud et al. (120021 ; l2005h . We will calculate the power spec- 
trum of this process directly from Eq. ( I A. 10b in the Appendix. 

We remark that for the amplitudes of individual events we 
assume the identical values (at each given radius). This re- 
striction is imposed only for the sake of definiteness of our 
examples; the formalism can deal with a distribution of am- 
plitudes. Indeed, we do not impose any serious constraint on 




Fig. 1. Illustrating the correspondence between the ignition mo- 
ments of the elementary events and the resulting light curve. 
The model is fully determined by a set of points in (t-r) plane, 
representing the pairs of ignition time t versus the time constant 
T of each event (panel a), and the elementary light curves with 
the profile I(t, r) - (t/rf exp (-f/r) 0(t) (panel b). The total 
light curve (panel c) is represented as a sum of the individual 
contributions. 



the model because the information about the level of the fluc- 
tuating signal can b e adjusted by setting the frequency of the 
events (lLehtoll989h . A simple demonstration of this concept is 
shown in Fig.[T] This plot illustrates how the model light curve 
arises from the elementary components. Naturally, we can ap- 
proach such decomposition from another angle, by investigat- 
ing how the total light curve can be expressed in terms of some 
basic profile. It is important to realise that, for the purposes 
of our present paper, light curves are of secondary importance. 
Instead, our calculations allow us to proceed from the distri- 
bution of the onset times and the characteristics of individual 
flares directly to the power spectral density, which stands as the 
primary characteristic of the source signal. 

Equation (|2]i represents a very general class of random pro- 
cesses. By applying the Fourier transform, we find 

2 sm(Ta)) v~' r 



rT[f(t)](co) 



xF(t-6pj,rj)](c^), (3) 

where • denotes the convolution operation. In Eq. (|2]), we sum 
together a set of all events (this infinite sum could in general 
pose problems with convergence, however, as we will see later, 
we can restrict the sum to a finite set of events without any loss 
of generality). The Fourier transform of a single event, I{t - 
5j, ^j)F{t - 6pj, rj), is then 



r[lit-6j,^j)F(t-6pj,rj)]ico) 

= e-''-'^r[l(t,^j)] * r[F(t + tpj, rj)](aj). 



(4) 



Periodical function F{t, r) can be now expanded in a series, 
F(t,r) - (2;7r)"' 2 Q(r)exp[/A:Q(r)f], and the expanded form 
substituted in the incomplete Fourier transform of /(f). 
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Knowing the incomplete Fourier transform of f(t), we can 
calculate its squared absolute value and perform the averaging 
over all realisations of the process. This can be simplified by 
assuming that every single event quickly decays. In principle, 
between -T and T the process is influenced by all events ig- 
nited during the whole interval (-00, T), however, because of 
the fast decay of the events, the interval can be restricted to 
{-{T + C), T), where C is a sufficiently large positive constant. 
In other words, every realisation of the process f{t) on {-T, T) 
can be described by a set of points in (k + 4)-dimensional space 
(tj, tij, tj,j, rj, ^j) with -(T + C)< tj < T. 

Values of the initial time delay and phase are functions of 
initial position of each spot, i.e. 



tdj = St(rj, <Pj), fp 



+ fd./. 



Q(0) "^- 

Fourier transform of the resulting signal is then 



(5) 



r[l(t - tij, ^j) Fit - taj + t^j, rj)\(cj) = Yj ^'t^'") ^"'* ^*' (6) 

k=-co 

where Tk = 'FUit - 5t{r, (p), ^)](w - kQ.(r)) is the Fourier trans- 
form of the event light curve, corrected for the time delay. 
Every realisation of this process is completely determined by 
the set of points [tj, (f>j, fj,^:). 
Defining the function 



s(t, (p, r, ^; w) 



we can write 



2 sin(rw) 



a> 



CO 



^=-00 



xr[I{t-6t,0]i<^-kilir))^ 



(7) 



2 



Due to Campbell's theorem ( lA.lSI l, 

AxA' 

X i*(f', 0', /, ^'; w) m2(f, (p, r, ^, r', 0', r',^') dA dA', (9) 

where m2 is density of the second-order moment measure M2(.) 
corresponding to the random point process of {tj, (f>j, rj, ^ ), A 
is a Cartesian product of four sets representing the domains of 
definitions, i.e. 



A = {-(T + C), T) X (0, 111) X <r„ 



K>xS. 



(10) 



Now we can perform the limit of T —> 00, as given by 
Eq. (lA.lOb in the Appendix. It can be shown that the result is 
independent of the value of C. Therefore, to obtain an explicit 
formula for the PSD we need only to specify the form of mii.) 
in Eq. (|9]l. Hereafter, we will show how to proceed with this 
task. 



2.3. Model driven by the Poisson process 

To start with a simple example, we assume mutually indepen- 
dent events with uniformly distributed ignition times. In other 
words, in this subsection we assume that there is no relation- 
ship among different spots - neither in their position nor in the 
time of birth (spots are statistically-independent). The intensity 
and the second-order measure are 



Mgi(df) = Adt, 



(11) 
(12) 



Mg2(dfdf') = [a^ + A6{t-t')\dtdt', 

where A is the mean rate of events. Other parameters are treated 
as independent marks with common distribution G{d(f>drd^). 
The second-order measure is 

M2(dfd0drd^df'd0'd^') = [a^ G{d4> dr d^) 
X G(dcf>' dr' df ) + AG(d(f>drd^) 
x6it-t')6{4>-4>')6(r-r')6(£-^')]dtdt'. (13) 

The result of integration (|9]l can be simplified for T 



The t ask reduces to evaluation of two limits (see iPechacek 
20081 for details). After somewhat lengthy calculations we ob- 
tain a general formula for the power spectrum: 

S(aj) = 47t\ Yj I ^k(r) c*,(r) e'('-*>^ Tk T,* G(d<^drd^). (14) 



kj=- 

Here, we remark that the general relativity effects are in- 
cluded in the Fourier coefficients Q(r). Knowing them in ad- 
vance (i.e., pre-calculating the sufficient number of the coeffi- 
cients that are needed to achieve the desired accuracy) helps us 
to produce the PSD efficiently. But we will start by neglecting 
these relativistic effects, so that we can clearly reveal the im- 
pact of the intrinsic timescales of individual spots and the form 
of the avalanche process. 



(8) 2.3.1 . Example 1 : Markov chain 



Let us consider box-shaped events with exponentially- 
distributed life-times, i.e.. 



I(t,T)^IoX(a,T)(t), ((T)^fie- 



(15) 



where x is the characteristic function (x - ^ for < f < t, 
X = otherwise); ^(t) is the probability density; and /q is a 
constant (we will set /q = 1 for simplicity). This again rep- 
resents a process of the type (|2|i, but at the same time it can 
be con sidered as a continu ous-time Markov chain with discrete 
states (ICox & Milleiill965h . The process PSD is then given by 
Eq. ([Ull with 

A' 



G(dTdr) 



-/"" 



drdr, F(r,r)=l. 



(16) 



In this example, by setting the transfer function F{t, r) = 1 
we completely "switch-off" the periodical modulation and the 
relativistic effects. 

Coefficients Q(r) are given by the relation 

27r/n(r) 

Q(r)= f F(f,r)e-*"W'df, (17) 
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which for the constant function F(t, r) leads to co(r) - **) 
(2n)-^ = 0.159, Q(r) ^ (k > 0). Fourier transform of the 
profile function is 

1 - e-'"^ 



ia> 

,2/ 



r[iit,T)](co) = 

\T\I(t,T)\{aj)\ = . 

Substituting into the general formula (fl4l l we find 

/■max oo 

4nAj\co(r)f J\r[l(t,T)]( 




S(aj) 



m 



dr dr : 



2A 



(20) 



Notice that the theory of Markov chains is usually formu- 
lated in_anotherwa}^_independent of our previous calculations 
(see lCox & Milleij[l965l) . One can thus obtain the PSD of the 
Markov process by two different approaches and verify the re- 
sults. Such comparison gives the same result, as expected. 

3. Introducing a relationship among spots 

The assumption that spots are statistically independent seems 
to be a reasonable (first) approximation. However, the actual 
ignition times and spot parameters should probably depend on 
the history of a real system. The statistical dependence among 
spots can be introduced in several ways. In this section we dis- 
cuss different models where an existing spot gives, with a cer- 
tain probability, birth to new spots. In this way a single spot at 
the beginning can trigger a whole avalanche of its offsprings. 
This avalanche can be in principle of arbitrary length, although, 
to obtain an infinitely long avalanche with non-diverging rate 
of new spots one would have to fine-tune the parameters. In or- 
der to avoid the unlikely fine-tuning and to obtain a stationary 
process, we will assume the occurrence of many spontaneous 
spots distributed by the Poisson process that keeps triggering 
new avalanches of finite duration. 

The first example can be called "Chinese process". By defi- 
nition, an existing spot gives birth to exactly one new spot with 
probability Q < if/ < 1 . In other words, every event produces 
at most one offspring. The spot of the kt\\ generation is always 
ignited later than the spot of the (k - l)th generation. As men- 
tioned above, spontaneous spots arise randomly, according to 
Poissonian process. In the simplest version of this model, de- 
lays between the parent spot and its lineal descendant are ran- 
dom values obeying the probability density pit). 

More generally, every spot can deliver n new spots, where 
« is a random value with Poisson distribution and the mean 
V. We describe this situation in terms of (i) standard Hawkes' 
process, and (ii) the pulse avalanche model. The temporal dis- 
tribution of new spots is no w governed eith er by the function 
p.{i) of the Hawkes' process (iHawkesI 1971 ), or u (t, r | fo, T(\) in 
the case of avalanches (Poutanen & F abianll 19991) . Spots of dif- 



ferent generations can appear at the same time. 

The difference between the Chinese process and the 
latter two processes is schematically sketched in Fig. |2] 



Fig. 2. Distinguishing between the Chinese process (a) and the 
Hawkes processes (b). In both panels, points represent ignition 
times of the spots. Each sequence starts with a spontaneously 
generated parent spot (open circles) and it continues with sub- 
sequent secondary ones (filled circles). Arrows symbolise the 
parent-daughter relation. The difference between the two sce- 
narios is described in Sects. l4n and 14.21 respectively. Within 
the schematic level of this graph, the pulse avalanche process 
(see Sect. 14. 3b belongs also to case (b). 



Mathematically, all three examples belong to the class of clus- 
ter processes. 



3. 1 . The cluster processes 

Point processes are characterised b y the generating functional , 
§[■], which is defined by its action jDalev & Vere-Jonesll2003h 



g[h(x)] = e[ ]^ h{xd 

x,esupp|A') 



(21) 



where h(x) : /Y — > C is a function satisfying the condition 
\h(x)\ < 1. 

The functional ^^[.1 satisfies various relations which can be 
derived in close analogy with the theory of generating func- 
tions of random variables. For our purposes it will be useful to 
expand 0[.] into a series in terms of factorial measures. 



g[l+T]] = l+Yjl j](xi) . . . T](xk) M[^(dxi . . . dxk). 



(22) 



The cluster processes consist of two point processes. The 
first one is connected with the counting measure NdA) and de- 
fines the central points y of the clusters. The second process 
spreads new points around the central point according to the 
random measure N(B\y). The complete counting measure is 
thenN(A)^ J^N(A\y)N,(dy). 

Let g [h(x) I y] be the generating functional of a cluster with 
the center at y, 



g[h(x)]^E^jg[h(x)\y]N,(dy) 



(23) 



If the cluster center process is Poissonian, the latter formula 
simplifies to 

@[h(x)] = exp { - J (l - ^ [h(x) I y] ) Ae(dy)}, (24) 
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where Ac(dy) is the intensity measure of the center process. 
One can expand the functional ( l24l i in terms of the factorial 
measure of the cluster, 

oo p 

^[l+77ly] = 1+^ I T](xi)...T](xk)Mit](dxi...dxk\y).(25) 

Putting Eqs. ( l22l i and dZST l into Eq. ( l24l i and collecting the terms 
with the same order of ri(x), we find 

M[i](A)=jM[i](A|3;)Ac(dy), (26) 

X 

M[2](A X B) = rM[2](A xB\y)A, (dy) + M[i](A) M[i](B).(27) 

A' 

For a stationary process, the intensity measure stays con- 
stant, Ac(d;t:) - A^dx, and gives the density of the centers. All 
factorial moments are shift-invariant in the sense 



mik](xu ..., Xk\y) = nimixi - y,. . . , Xk-y\0), 



(28) 



where m[y(. | y) is density of M[y(. | y). As a result of the shift 
invariance, density of the first-order moment nii must be also 
constant. Moreover, we can always choose y in Eq. (|28T l equal 
to one of Xi, so the function mj^j depends on only k- 1 variables. 
We define the reduced factorial moments. 



mik](uu- ■ ■ , Uk-i) = m[i:](0, mi, . . . , Uk-\), 
and from Eq. ( |27] i it follows that 



m[2](M) = /i I m[2](3',3'H-M|0)dy-l-m 



(29) 



(30) 



Stationarity of the process implies that the second-order mea 
sure density depends only on the diff'erence of its arguments, 

m\2\(t, f) - OT[2](f - f) - C(t - f) + TO), 



where c(f) is an even function. 

The second-order measure of a marked cluster process is 

M2(dtd^drd^dt' d0'd^') = [(m^ + c(t - t'))G{d^drd^) 
X G(d0' dr' d^') + my G{d(pdrd^) 
x5(t-t')6{cl>-<p')5{r-r')6{^-0]dtdt', (32) 

almost identical with that of the Poissonian process. There is 
only one additional term associated with the function c{t). The 
resulting spectrum is given by a somewhat lengthy, but still ex- 
plicit formula. We find stationary cluster processes to be par- 
ticularly promising as a general scheme, encompassing a broad 
range of models as special cases. The PSD is 

S{aj) = S lio)) + An^mi 

OO p 

X 2 J e'('-*'%(r)c;(r)5",(w)5";(w)G(d<^drd^) (33) 

k,l=-oo ^ 

with 5 i(w) = 47r25p(w) Q,^{(jj) Q.*'{cj) and 
Q.,{c) = £ re-'*^Q(r);ri(w)G(d0drd^). 



(34) 



The reduced quadratic factorial moment appears in the for- 
mulae for power spectra of cluster processes expressed in terms 
of Sp((jS) = T[c(t)]{co) = T[m[2](t)-mlj{co). This is di- 
rectly related to the two-dimensional Fourier transform of the 
quadratic factorial measure as 



Sp{co) - A m[2](tj, -co I 0), 
where 

m[2](w, w' I 0) = re'*^'"+-^''-"'>mp](xi, X2 \ 0) djcj djC2. 



(35) 



(36) 



We will discuss the form of Spico), the expression for mi, and 
the resulting PSD in diff'erent situations. But before that, we 
still need to show how the model properties can be detailed in 
term of marks. 

3.2. Marks as a way to specify the model properties 

Until now the variability patterns have been restricted only by 
very general properties of the assumed process. This means that 
the model is kept in a very general form. However, formulae 
(fT4l l and ( l33T l are too general for any practical use. Their main 
value is after defining special cases. Then these formulae can be 
readily applied to derive the analytical form of the PSD. Such 
special cases are conveniently defined by means of marks. We 
discuss possible choices of the mark distribution, G(dr d0 d^). 
We can simplify the situation by assuming axial symmetry. 
Therefore, all statistical properties should depend only on the 
radius (the azimuthal part of G is constant). The distribution of 
marks has now the form 

G(dr d0 d^) = 1 ^r(^ I r) p(r) dr dcf>d^, (37) 

27r 

where (j^{^ \ r) is the reduced probability density of the remain- 
ing parameters. To illustrate this case more c learly, we will now 



(31) examine th e phenomeno logical model of lAbramowicz et al 



(1199 ih and IZhang & Bad (Il99lh . which can be considered as 
representation of the spotted accretion disc. 

3.3. Example 2: a spotted disc 

Let n(r) be an average number of spots at radius r, each of 
them shining with the average intensity /m(0 for average dura- 
tion Tuir). Let us further assume that these characteristics scale 
with the radius as power laws: 



n(r) = A„ r"" , TMir) = A, r"^ , Mr) = Aj r" 



(38) 



(here, as and As are constants). This setup falls within the cat- 
egory of models described by Eq. (|2]l, with the profile function 
I(t, r, t) and the mark distribution G(dr d(p dr). The first two re- 
lations (|38] | prescribe the conditional marginals of G, i.e.. 



xW 



n(r) = Ap(r), Tuir) 



= f r^R(T|r)dT, 

T„,i„('-) 



(39) 



where the integration limits can depend explicitly on r. The 
third equation determines the amplitude of the profile function. 



k=-oo' 



I(t,r,T)^A,r-"g(t,T)eit). 



(40) 
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we rewrite the PSD formula ( fT4l i in the final form 



The dependency on t is not determined uniquely. In order to 
o btain an exphcit form of G we have to go beyond the model 
of lAbramowicz et alj (Il99lb by assuming the distribution of r. 



Sioj) = 4n 



"Mdk-ir) 



(41) 



(42) 



(43) 



The mean, Tyiif), must satisfy Eqs. ( [38] ) and ( |43] ). The choice 
of 



UT\r)^K(r)T-''. 




The normalisation constant then equals 




F(^^ - ^~P 








and for tm(0 we find 




Tuir) = Kir) t' ''dr = — 





n/j I '^k('-)c'',( 
X JTkico) r;(aj) ^r(^ I r) d^p(r) dr. 
We find the coefficients by direct evaluation, 

2n 



(54) 



d„(r) 



1 fe'" 

27: J 



1 + Q.(r) 



d6ti^(y),r) 



d(p 



(55) 



T ■ (r) = C ■ r"' T Cr") = C r"' 



(44) 



is consistent with both equations. Constants Cmi„, Cmax and A^ 
are bound by the relation 

Ar(2 - p){C]^^ - C^;^) = (1 - p)(ci7x - C^D (45) 

Because Cmin and Cmax are positive, we can write 

Cmin = C, Cmax = jC, (46) 



C ^ Ar 



K(r) 



2-p y^-P - 1 
1 - p y^-P - 1 ' 

■i^^(Cr'^0'^-V 



yi-P - 1 

Therefore, by substituting back to Eq. dTO . we obtain 

^R(T\r)p{r) - 



(47) 
(48) 



with y - (f) + Q(r) 5f(0, r). We note that the term ( fSSl ) corre- 
spon ds to the effect of delay amplification in terminology of 
Dov ciak et alj ( 120081) . It influences the observed signal from a 
source moving (i.e., orbiting) near a black hole. 

4. Results for the model PSD 

4.1. Model driven by the Chinese process 

Let us denote if/ the probability that an existing spot generates 
a new one, and qt the probability that a family of spots con- 
sists of exactly k members. The value q^ obeys the geometrical 
distribution, q/t - i/'*(l - iff). 

We interpret probability density p{t) of the delay between 
successive spots as a mean number of first-generation spots that 
occur at the ignition time f > 0, where f = is a moment of 
ignition of the parent spot. Analogically, p{t) • p(t) is the mean 
number of second-generation spots. For a sequence of k spots, 
we obtain the intensity measure 



\ ^^^ mm / 



- 1 



mi(t\0,k)^J]p*\t), 



(56) 



7 _ i-p _ 1 \'''^' 

X 1^ — ^^„ -A J T--p/p-i)«.+<^„^ (49) 



7=0 



,1-/7 y^-P - 1 
where the definition domain is a set 

E = {(r, ^r)) I r e <rmi„, rmax), ^r) e C/''<0, y}} . (50) 

We can set p - 1 as a specific example. This value of the 
power-law index is special in the sense that neither short nor 
long timescales dominate the PSD, as we see from Eq. (l42l i: 



where p*''{t) is the ^th convolutionary power of p{t). We can 
write mi in the form 



miit\0)^Yjqkmi(t\Q,k). 



(57) 



k=0 



Convolution of two functions can be calculated via the Fourier 
image. Defining p{cj) = T [p(t)] (to) we get 



Iny 1 

C^Ar^, Kir)^-—. 
y - I Iny 

By substituting back to Eq. ( |4TI ) we obtain 
^r(t I r)p(r) = 



(51) m 



i(wio) = (i-^)y/y/7V) = ^ — -T--- (58) 



From this we find 



1 / Q-n + l ffn + U 

lnr(^max -'"mm j 



r"-'T-\ 



(52) 



Knowing the concrete form of the distribution of marks, 
we perform the integration over ^. As l(t, r) does not explic- 
itly depend on azimuthal angle, the integration is simplified. 
Denoting 



OO 

/ 



CO 

I 



mi(f|0)df = mi(w|0)|„=() = 



1-r 



tmi(t\0)dt = -/--— mi(w|0)|a;=o = 
dw 



«A 



ii-ifry 



(59) 



■E[t]. (60) 



In 



dn{r) = ^ ' ^' 



27r J 



m[<t,+n(r)6t{r,il>)] 



d(/>. 



(53) 



The meaning of integral ( |59] ) is the average number of spot off- 
springs in the whole avalanche. The meaning of the last integral 
is the average duration of the avalanche. 
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Calculation of the quadratic measure is a less intuitive pro- 
cedure. We start from the generating functional (I^TT i of the pro- 
cess. 



For a stationary process, the first-order moment density is 
constant. Averaging the relation ( |66] l we find 



^[/!(0|0] = £?//z(0)^,[/!(0], 



mi 



(61) 







(68) 



/=() 



where 0i[.] denotes a generating functio nal of finite renewal 
proce ss with I points (see chapt. 5 in iDalev & Vere- Jones 
' 2003h . 

Substituting h{t) -I- rjit) in the expansion (l22l i. we obtain 
the Fourier image 

i/r[/7(w) + p{u') - 2l/r/7(w)p(w')] 

OT[2](W, w' 1 0) = r.(62) 

[l - ilfp{cS)\[\ - ilfP{u')\[\ - ijfPiu + w')J 

This equation allows us to find the second-order measure. The 
PSD is then given by Eqs. (O-dHi with 



The meaning of v is the mean number of the offsprings. Clearly, 
it satisfies the normalisation ^ ^{x)Ax = v < 1. 

The generating functional of the cluster of the Hawkes pro- 
cess fulfils the integral equation, 

g[h{x) I O] = h{0) exp { - J (1 - ^ [h{x) I y]) fx(y) dy}. (69) 

— oo 

Substituting h(x) = 1 -t- rjix) and expanding both sides into the 
series (IZST i we find 



nil — 



l-i/f' 



Sp(cj) 



2Aifr[ne[p(aj)] + ^\p{aj)\^] 



(l-^)[l-2if,%e{p(aj)] + i/^^p(cj)\2]' 
Equation ( l64b can be cast in the form 

|m[i](w|0)|2(l-^2|p(w)|2)-l 



(63) 
(64) 



OO 

m[i](f|0) = I mii]{t\y)fi{y)dy + 6{t), 

— OO 
CO 

m[2] (t, f' I 0) = I mp] (t, f I y) p{y) dy 



(70) 



+m,u(t\0)mm(t'\0)-6(t)6(t'). 



(71) 



Sp(a>) = A ■ 



1-«A 



(65) 



To complete the calculation we solve the integral equation 
jTOl l for m[i](x 1 0). Because this is a linear convolutional inte- 
gral equation, it can be solved efficiently by using the Fourier 
transform; 



4.2. Model driven by the Hawkes process 

Hawkes' process (lHawkeslll971tlBremaud & Massouliell2002l) 
is more complicated than the previous example because it con- 
sists of two types of events. First, new spots are generated by 
the Poisson process operating with intensity A (let us denote fo 
the moment of ignition). Second, an existing spot can give birth 
to new one at a later time, t, according to the Poisson process 
with varying intensity jj.(t - f())lj 
The mean number of events is 



OO 

/ 



//(«) = e-"^'//(Odf, 



OO 

m[i](w|0) = re-'"'m[i](f|0)df=y-4 



liiui) 



(72) 



(73) 



For the Fourier transform of the quadratic factorial measure we 
find 



m[2](w,w' |0) - 



OT[i](w|0)m[i](<y'|0)- 1 



OO 

1(0 ^ A+\^yi{t-t^^ A+ I /i(f - X) N(dx). 



(66) 



In analogy with Eq. ( I60b we define the characteristic time 
of avalanche fa- It can be proven that fa is related to the charac- 
teristic time of the infectivity ti. 



1 - jj{<jj + Cl)') 

Again, the PSD is given by Eqs. 0311, OU, and ^8^ with 

\nnii(co\Ot - I 



Sp{cij) = A ■ 



1-v 



(74) 



(75) 



j'fmi(f|0)d? ftij(t)dt 

_ ^ _ V 

"^ ~ 1-v^ ~ 1-v 

/TOi(f|0)df /MOdf 





fi. 



(67) 



Comparing this equation describing the Hawkes process PSD 
with the corresponding Eq. ( |65] l for the case of the Chinese 
process, we reveal a subtle difference between the two mecha- 
nisms. It turns out that the high-frequency limit is identical for 
both of them, however, the difference grows as one proceeds 
towards the low-frequency end of the PSD domain. 
For the exponential form of infectivity measure. 



' Mathematically, a very similar model has been employed 
to describe the propaga tion of diseases through a population 
( IDalev & Vere-Jonesll2003r) . In this context, the function ii(t) is called 
"infectivity", for obvious reasons. We will adopt the same name for it, 
although the medical connotation is irrelevant here. 



iu(t)^vo-e-°''0(t), 

we obtain the expUcit expression of 

m[i](w|0) = 1 + 

Sp(CL)) = 



cr( 1 - v) -I- ia> 
v(2-v)o-2 A 



(1 -vVo-^ + co^ 1-v 



(76) 

(77) 
(78) 
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Fig. 3. Power spectra from the spot model in which the birth and duration of spots are governed by the market Hawkes process 
with the exponential infectivity ( |76] |. Values of a are given on top of the plots. Frequency is given in geometrical units (see the 
text for its scaling to physical units). In each frame, the four curves correspond to different values of v; from bottom to top v = 0. 1 , 
0.5, 0.9, and 0.99. The mean number of spontaneous spots has been set to A = (InyK The relative normalisation of these curves 
scales as the mean number of all spots, i.e. proportionally to 1/(1 - v). Left: the case when all profiles I(t) are identical; t = 1, 
i.e. ^(r) = 6{t - 1 ). Right: the case when the life-times of spots are distributed uniformly, i.e. according to ^(r) = 1 /(T,nax - T"min), 
between Tmin = 0.01 and rmax - 1- 



where o" > is a constant. 

Figure [3] shows the resulting PSD of this model in a loga- 
rithmic plot of ojS (co) versus a>. Here we can study the occur- 
rence of break frequency where the PSD slope changes depend- 
ing on the model parameters. Light curve profiles of individ- 
ual spots were chosen as exponentials, /(f) = /q exp(-f/T) 0{t), 
where r is random value with probability density ^(r) and the 
mean f. The characteristic time of infectivity, t, - cr"', is set to 
be ti - a f. In general, we can identify two characteristic fre- 
quencies in the PSD. The first one corresponds to the charac- 
teristic frequency of the profile I{t) (in our case this frequency 



is given by 1 /f), the second one is given by the characteristic 
frequency of the avalanches (1/4). 

We remind the reader that this plot (as well as the subse- 
quent figures |4]-(6]l does not include general relativity effects; 
they will be recovered later in the paper. This is merely to 
simplify calculations - the relativistic effects complicate the 
derivation of the analytical formula for the PSD and it is some- 
what difficult to distinguish them from the intrinsic properties 
of the signal. 

The typical form of the model PSD resembles a Lorentzian 
or doubly-Lorentzian profile at central frequencies. In most 
cases, there is only one break in the spectra corresponding 
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Fig. 4. Graphs of the model PSD as in the previous figure, but for the case of Hawkes process with the power law infectivity 
(l80t . Parameter values fi are given on top of each panel; values of v are as in Fig[3] All spots have the same (exponential) profile, 
as in the previous figure. Fory6 > 1.5 we notice the new profile resembles a power law that develops in between the two peaks. 
However, this situation is rather rare within the parameter space of our models. By parameter tuning - e.g., by setting v — > 1, 
which means enhancing the contribution of avalanches while suppressing the importance of spontaneous parent spots - the extent 
of the flat part of the PSD profile (like the one seen for^S = 1 .5) can be stretched farther towards small frequencies. 



to the characteristic frequency of the profiles I{t). When two 
peaks are visible, as we notice from Fig. [3] and Eq. (fTSl l. two 
characteristic frequencies occur The upper frequency is scaled 
to unity in our figures; this peak corresponds to the e-folding 
time of the exponentially decaying spots, i.e. I{t) oc exp(-f/T)|3 
The lower frequency peak is then set by a combination of two 
timescales. 



1 -V 



a 



(79) 



^ A single exponentially-decaying pulse is described by function 
I{t,^) = In e^'l^6(t). Notice that in this simple case the only parameter, 
T = ^, has a meaning of the decay timescale. 



where f s E[t]. Notice that the model is sufficiently complex, 
and so one cannot easily disentangle the two timescales in any 
straightforward way directly from the PSD. 

Another natural choice of infectivity is a power-law func- 
tion of the form 

^« = (7^^«' ^''^ 

where fi, K, and L are positive constants satisfying 

/3>l, K^v(J3-l)LP-K (81) 

One can thus ask how the PSD form depends on the assump- 
tions about the form of the infectivity function //. For every 
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a= 10, v = 0.9 





10"^ 10"^ 10"^ lO"'^ 10"^ 10"^ 10"^ 10° 10^ 10^ 10^ 



Fig. 5. Decomposition of the total PSD curve, St(co), into a product of two terms which are responsible for the two peaks in the 
final profile (see Eq. (l83Tl). Left: this case corresponds to or = 10, v = 0.5 curve shown in the bottom-left panel of Fig. [3] Right: 
this case corresponds to /? = 1 .9 and v - 0.99 curve in the middle-right panel of Fig. |4] Similar behaviour of the PSD profile is 
typical in our model, although in some cases only one peak dominates the spectrum whereas the other one is suppressed. 



b > the function ( l80t satisfies /u(Olz,=i = fe"' yu(f//7)|^^[. 
Therefore, we can set L = 1 without any loss of generality 
(the value of L can be recovered by a simple rescaling of the 
result). 

We note that the characteristic times t-^ and fi, as defined 
in Eq. (l67b . diverge for the infectivity (l80b and 1 < /? < 2. In 
this case, fa does not exist. Again, we derive an analytical form 
of the PSD for the adopted infectivity function. The procedure 
is entirely analogical as above, but we do not give the explicit 
form of Sp{a)) because the final formula is rather complicated. 
The resulting PSD curves with the power-law infectivity are 
plotted in Fig. m 

It is also interesting to note at this point, how the peaks 
of the PSD move when the model parameters are shifted. For 
example, see the panel o- = 1 of Fig. [3] Although in this case 
the two timescales are equal to each other, changing the other 
parameter, v, from 0.1 to 0.99 brings the peak over two orders 
of magnitude. In other words, the maximum of the PSD can 
appear at a frequency lower than the inverse of the spot decay 
time. The frequency shift of the peaks is again given by factor 
fa/f, as explained in Eq. ( |79] |. 

In order to understand better the behaviour of the PSD 
peaks, we rewrite Eq. ( l33l l in the form 



S(cj) = mi E[|5"[/](o.)f] + SAco) \E[r[I](co)f 



(82) 



(i.e., relativistic effects neglected). Equation ( |82] | can be simpli- 
fied by assuming that all spots have the identical, exponentially 
decaying light curves, I(t, t) - /o e^'^'^ 9{t). We obtain 



Siu) - S\^{Lj)(m\ -i-5p(w)j, 



(83) 



where Si^ico) - t 1(1 + w r ) is the Lorentzian PSD, and m\ 
is a constant (see the discussion following Eq. (l28ll). In this 
simple case the two terms, SAcS) and 5p((jj), give rise to two 
comparably significant peaks, although their amplitudes are 
generally diff'erent. The situation is shown in Fig. |5] The left 



panel is particularly transparent because it shows the Hawkes 
process with exponential infectivity, for which Sp{u) adopts 
the Lorentzian shape, dominating the spectrum at frequencies 
where Si^icS) ^ 1. Therefore, Eq. ( l83T l is eff'ectively a sum of 
two identical Lorentzians. However, this example does not ap- 
ply to a general case when the amplitudes of the PSD peaks 
can be very different. In the right panel, we note that S p(a») it- 
self dominates the whole spectrum and has a complicated (non- 
Lorentzian) shape. 

4.3. The pulse avalanche model 

The pulse avalanche model was discussed in the context of var- 
ious astronomical objects whose light curves exhibit signs of 
stochasti c behaviour. They are, namely, gamma-ray bursting 
sources dStern & Svenssonlll996h . As a framework to describe 
the timing cha racteristics of accret i ng bla ck holes, the model 
was studied by IPoutanen & FabianI (119991) . Details of the pro- 
cess are diff'erent in those two papers and our description is 
closer to the latter one. 

The basic properties of the pulse avalanche model can be 
summarised as follows, (i) The observed signal consists of 
pulses of the form I(t, r) - Iq g{t, r) d(t), where r is their char- 
acteristic duration, (ii) Each pulse gives birth to b baby pulses; 
the number of baby pulses varies, obeying the Poisson dis- 
tribution P(b) - v^ l{b\e^) with the mean value v. (iii) The 
baby pulses are delayed with respect to the parent ones by Af, 
which is a random value with exponential distribution, P(Af) - 
(Q'To)"'e"*'^''"^°'. (iv) Some pulses occur spontaneously, ac- 
cording to the Poisson process operating at the mean rate A. 
Finally, (v) temporal constants, r, for all pulses are mutually in- 
dependent and drawn from the same distribution function, ^(t). 
Such a process is clearly of the form (|2]i. 

The underlying point process is a cluster process operating 
on the set Rx(Tniin, Tmax) with properties similar to the Hawkes 
process. (In fact, we will see below that the Hawkes process 



12 



T. Pechacek, V. Karas, & B. Czemy: Hot-spot model for accretion disc variability 




1 
1 
1 
1 
1 

1 ^ 

2" 1 

1 
1 
1 
1 
1 



10-' 




10'^ 10"^ 10"' 10° 10^ 10^ 10^ 





Fig. 6. The pulse avalanche model for the uniform probability density ^(t), constrained by limits rmin = 0.01 and r^ax - 1- 
Left: the PSD calculated from the formulae (l33b and ( l86b . Right: the ratio between the PSD of the pulse avalanche process to 
the Hawkes process; in the latter the exponential infectivity /i(f) is related to the corresponding infectivity of the pulse avalanche 
process by ju(f) - fi(t \ f). 



can be considered as a special case of the pulse avalanche pro- 
cess.) The center process is the Poissonian one, and its intensity 
is /lc(f, t) = /1^(t"). The clusters are driv en by the Poissonian 
branching process JCox & Millen Il965h with the parameter 
measure 

f - f 1 



//(f, T I fo. To) = exp 



I- 



9(1 - to). 



One can prove, by direct integration, that 

E[Af I fo, to] = - I I (f - to) fi{t, T I to, to) dr df = aro, 

—CO Tmin 
CO r^ax 

e[/7 I fo, To] = I I nit, T I fo. To) dr df = V. 



(84) 



(85) 



In analogy to Eqs. ( iTOl i and ( fTTI ) one can derive a set of linear 
integral equations for the first and second order moments of the 
cluster process. The chosen form ( l84b of the kernel fi{t, T\to, To) 
greatly simplifies the task, as it can be again transformed into a 
convolutory form and solved in the Fourier domain. 



The resulting functions ihm and m[2] are complicated and 
their inverse Fourier transforms can be directly found only in 
very special cases, e.g. by assuming ^(t) - 5{t - a). In that 
case the pulse avalanche model is transformed into the Hawkes 
process with the exponential infectivity measure. 
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Fig. 7. Upper panels: Power spectra from tlie Poisson-driven spot model, calculated for a thin accretion disc extending between 
radii r = 6M and r = lOOM (in geometrical units), for three inclinations Oq. The purpose of this plot is to demonstrate a general 
agreement between the PSD calculated from the model light curve and from the analytical formula. The wiggly (red) curve is a 
result of direct numerical simulation, including the relativistic effects. The smooth (blue) curve is the analytical result, derived 
from Eq. (fT4l i by specifying the probability density function ({t) oc 1/t. The vertical (magenta) lines represent the Keplerian 
orbital frequencies, 0(r), at the inner and the outer edges of the disc. Lower panels: Power spectra from the Hawkes process with 
the exponential infectivity (a = 7, v = 0.8). The analytical curve was calculated by using formula (1331) and assuming the same 
probability density function ^(r) oc 1 /r. In all panels we set Tmin = 300, Tmax = 5000. 



Fortunately, we do not need the explicit form of the mo- 
ments to find the final formula for the PSD. This result is pro- 
vided by Eq. ( l33b with 

5i(w) = 4n^ f fspioj, T, t') Q^(oj, t) Q*(aj, r') dr dr' (86) 



general form of the PSD profiles exhibits very similar trends. 
We see only quite minor differences in the limiting slopes of 
the PSD and in the overall PSD shape in the middle range of 
frequencies. Th erefore, we do not show these plots here (see 
|Pechacekll2008l) . 



and 



Sp(u,t,t') = ^(T)f(r')f/i(w)|/2(w)l' 

1 - V L 

+^ — /2(w) + , " , /;(«)], (87) 

1 — laTOj 1 + lar CO J 



/i(c^) 



Jr„„„ 1 + a 



(iy) 



7 9 7 



dy, 



fiioj) = 



Jr„„ 1 



Ay . 



+ 1 a toy 



(88) 
(89) 



We show several graphs of the resulting model PSD curves 
in Fig. |6] The overall form of the PSD resembles the previous 
examples derived for the Hawkes process, although they are 
not the same. The two models - the pulse avalanche process 
and the Hawkes process - are identical in their high-frequency 
and low-frequency limits, but they differ from each other in the 
middle range of frequency. In order to demonstrate the differ- 
ence in a clear way we plot the ratio of the models in the right 
column. Fig. |6] assumes ((t) = l/(Tn,ax - T"min)- We also stud- 
ied the case of ^(r) = [Tln(Tniax/T"min)] ' and we found that the 



4.4. PSD with relativistic effects 

Relativistic effects influence the photon energy in the course 
of light propagation through the curved spacetime, towards a 
distant observer. In this paper we consider only the energy- 
integrated light curves, but even those are affected because the 
energy shifts modify the observed flux level. Furthermore, light 
rays are bent as they pass near the black hole, causing the light- 
focusing effect. In consequence the observed light curves dif- 
fer from their intrinsic profiles produced at the point of emis- 
sion, and this further complicates the decomposition of the PSD 
spectrum into elementary components. Therefore, in general 
the power spectrum cannot be rigorously expressed as a combi- 
nation of Lorentzians. However, in most circumstances the rel- 
ativistic effects are not very prominent - only a small fraction 
of rays passing very close to the black hole horizon and those 
crossing the caustics are affected. One expects that they cannot 
be ignored if the accretion disc extends down to the innermost 
stable orbit or if some non-negligible emission arises below 
that orbit. Also, high-inclination objects are affected more be- 
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Fig. 8. The analytical curves of the PSD profiles are plotted, with relativistic efi^ects included. In order to allow comparison with 
previous figures, we selected appropriate combinations of the model parameters: the infectivity a, the mean number of secondary 
spots V, and the spot distribution are the same as in Fig. [3] The inclination 0o and the spot distribution as in Fig. [7] The vertical 
lines again indicate the range of orbital frequencies corresponding to the assumed range of radii, 6M < r < lOOM, where the 
spots are distributed. The analytical form is quick to evaluate, hence it is convenient to obtain the PSD form for variety of different 
situations. 



cause in those cases the disc is seen edge-on and the intrinsic 
fluctuations of the emission are considerably amplified. 

Exemplary power spectra, including the relativistic effects 
are modelled in Fig. [T] where the numerical simulation is com- 
pared with the analytical result. The upper panels assume that 
orbiting spots are generated by Poissonian process; the lower 
ones show the PSD derived from Eq. (l33b . The assumed radial 
distribution of the spots was p{r) - pq{\- y/EJr) r"^, where po 
is normalisation constant. It turns out that the relativistic effects 
influence the final PSD especially at high frequencies and high 
inclinations. 

We remind the reader that frequencies in these plots are 
given in geometrical units (in physical units frequencies scale 
inversely with the mass of the black hole). We notice that a 
high-frequency part of the spectrum decays as a>S(a>) oc w"', 
whereas the break occurs towards lower frequencies. These 
plots provide us with graphical comparisons between the an- 
alytical form and the corresponding results of numerical sim- 
ulations. We note that the adop ted approximation of relativis- 
tic effects jPechacek et al.l2005h holds for moderate inclination 



(< 70 deg). It loses accuracy when the view angle becomes 
almost edge-on, although the main trend of the PSD remains 
unchanged. 

The main advantage of the analytical method is, obviously, 
in the possibility of obtaining a general form of the PSD, in- 
cluding the relativistic effects. We are ale to search systemati- 
cally through the vast parameter space of different models for 
which the model PSD can be explored across a wide range of 
frequencies. We take the advantage of this approach and plot 
variety of profiles in Fig. [8] Here, we omit the numerically- 
simulated curves of the previous figure, so the entire graph is 
constructed very efficiently. We assumed that spots are gener- 
ated by the Hawkes mechanism. 

Basic trends of the PSD shape are readily recognised. In 
particular, the curves have either one or two maxima, promi- 
nence and position of which changes with parameters. Notice, 
for example, the bottom right panel in which the PSD is almost 
flat over several decades of frequency, well below the typical 
orbital frequency in the disc (vertical lines denote the Keplerian 
frequencies at the edges of the assumed spot distribution). The 
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Fig. 9. The ratio R^{(jS) of PSD computed with and without relativistic effects. The difference is most visible around intermediate 
frequencies; Keplerian orbital frequencies at the outer/inner edges of the spot radial distribution are again indicated by vertical 
lines. The arrangement of panels and the assumed parameter values are the same as in the previous Fig. |8] In particular, the 
vertical lines indicate the minimum and the maximum orbital frequencies corresponding to the outer and the inner edges of the 
spot distribution. 



flattish part of the spectrum can be extended further, to lower 
frequencies, by enlarging the v parameter towards unity, i.e., 
by protracting the sequence of avalanches. Next, for large incli- 
nations we notice that relativistic effects produce a prominent 
bump. This feature occurs near the orbital frequency of the in- 
ner disc. Relativistic effects are the main cause of differences 
between this figure and Fig. |3] in which those effects were ne- 
glected. 

As already mentioned, general relativity affects mainly the 
high-frequency domain of the PSD, around the orbital fre- 
quency of the inner disc, where it adds power to the observed 
PSD. In physical units the relevant frequency generated at ra- 
dius r is Tll^ K 10"^(r/10/?g)(M/10^Mo) [Hz]. On the other 
hand, it does not influence the middle part of the spectrum, 
i.e. at frequency < ^".^^(rout), neither it changes the asymptot- 
ical form at far ends of the frequency range (where the PSD 
decays as power law). It has been argued that the additional 
signal is actually not prese nt in the data of MCG-6-30-15 
(IZycki & Niedzwieckil2005') . although the light curve from the 
long observation should reveal some excess. However, the sit- 
uation here is more complex because of the avalanches con- 



tributing power to lower frequencies. In fact, our model pre- 
dicts rather weak enhancement near the inner edge orbital fre- 
quency - unless the inclination is almost edge-on (which is un- 
likely). More power is therefore typically expected at moderate 
frequencies, lower than the inner-edge orbital frequency. 

In order to see at which frequencies the relativistic ef- 
fects are most important, we construct the ratio R^{(jS) of the 
PSD calculated with (denoted by superscript "gr") and without 
("cl") these effects taken into account: 



«s(w) = A^o 



5|Xw) 
5^'(w) ' 



(90) 



where the normalisation factor is No = Sp{0)/Sp{0). (Doppler 
and lensing effects are neglected in the calculation of Sp{aj).) 
Figure |9] shows that the PSD graphs are affected only by a 
constant shift at their low-frequency as well as at the high- 
frequency limits, i.e., the ratio /?s(<^) reaches a constant value 
at both ends of the frequency range. This is quite easy to un- 
derstand as at the ends ends of the frequency domain the rela- 
tivistic PSD acquires the same slope as the non-relativistic one. 
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Relativity shapes the PSD especially within the range of orbital 
frequencies of the assumed spot distribution. 



5. Conclusions 

We adopted the viewpoint that the variability pattern is deter- 
mined by the interplay among the bulk orbital motion, relativis- 
tic effects, and the intrinsic changes of the inner accretion disc. 
We concentrated our attention solely on the PSD characteris- 
tics. The spots have a certain kind of memory in our model. 

We gave several examples in which the PSD changes the 
slope and certain break frequencies. The frequency of the break 
depends on the interplay of model properties, i.e., the intrin- 
sic form of the spot light curves, which determine the individ- 
ual contributions to the total signal together with the avalanche 
mechanism. The location of spots on the disc and the inclina- 
tion of the source define the importance of relativistic effects. 

In some cases, a double break occurs and the overall PSD 
profile is then approximated by a broken power law. This is a 
promising feature in view of applications to real sources with 
accreting black holes. The broken power-law profile either re- 
sembles a combination of the Lorentzians or, in some cases, an 
intermediate power law develops and connects the two peaks 
across the middle frequencies. The change of the PSD slope 
is clearly visible and well-defined in some cases, though under 
most circumstances it appears rather blurry. The low-frequency 
limit of the PSD slope is a constant; the-high frequency be- 
haviour depends mainly on the shape of the spot emission pro- 
file, including the general relativity effects. In our calculations 
the emissivity was decaying exponentially and the slope of the 
PSD was equal to -2 at high frequencies. In between those 
two limits the intrinsic PSD is influenced by both the individual 
light curve profile and the underlying process. 

It is interesting to notice that the doubly -broken power law 
occurs only for certain assumptions about the intrinsic light 
curves of the individual spots or avalanches - their onset and 
the decay; in other cases the break frequencies are not well de- 
fined, or the broken power-law PSD is not preferred at all. We 
stress that if two Lorentzian seem to dominate the PSD (i.e., 
two peaks show up), it still does not necessarily mean that two 
single oscillation mechanisms operate simultaneously. Instead, 
it may well be the manifestation of the avalanche mechanism. 

We employed a general statistical approach to the variabil- 
ity of a black hole accretion disc with orbiting spots that contin- 
uously arise and decay. The origin and evolution of spots were 
described by Poissonian and Hawkes' processes, the latter one 
representing a category of avalanche models. We derived ana- 
lytical formulae for the PSD, Eqs. (fl4l i and ( l33T l, and we dis- 
cussed their limitations and accuracy. The main advantage of 
the analytical form is the insight into the properties and the 
fast evaluation that captures the main trends of the PSD shape. 

It is worth noting that the PSD does not maintain all infor- 
mation ab out the light cu rves that can be studied by Fourier 

Extensions have been discussed 



methods (Vioetal. 



1992 ) 



Our approach allows us to investigate the influence of the 
assumed mechanism, which describes the creation of parent 
spots and of subsequent cascades of daughter spots. In particu- 
lar, we can discuss the PSD slope at different frequency ranges 
and locate the break frequencies depending on the model pa- 
rameters. The relationship between the mathematical nature of 
the process and the PSD of the resulting signal is an inter- 
esting consequence of this investigation, as it provides a way 
to grasp and rigorously constrain the physical models of the 
source. Therefore we believe that the method that we described 
is very helpful for identifying the underlying mechanisms that 
shape the PSD in black hole accreting sources. 
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Appendix A: Spots as a random process 

A.1. Assumptions, definitions and preliminaries 

In this Appendix, we briefly introduce the formal mathematical 
approach and notation used throughout the paper We employ 
t he concept of random values on proba bility space, (Q, 2, P) 
(Kolmogo rovlll950l:ICox & MUleJll965h . where fi is the sam- 
ple space (i.e., the set of all possible outcomes of an experi- 
ment), E denotes the cr-algebra on fi, and P is non-negative, 
cr-additive measure satisfying conditions P(0) - 0, P(Q) = 1 . 
It is usually assumed on physical grounds that real signals sat- 
isfy all mathematical prerequisites. 

A real random value is a map X from D. to real numbers, 
X : O — > R. The distribution function F{x) and the probability 
density function f{x) are then defined by 



F(x) = P ({nr € n I X(m) < x}) , 



fix) = —Fix). 
Ax 



The mean value operator E acts as 
E[giX)\ = ^ gix)FiAx) = ^ gix)fix)Ax. 



(A.l) 



(A.2) 



Random values are characterised by their moments, defined 
as Mk = E[X*]. The above-given relations can be generalised 
to higher dimensions by introducing a A:-dimensional vectorial 
random value, X - iX\ , . . . , Xk), the distribution function Fix), 
and the probability distribution 



f(x) 



d 



d 



dx\ dxk 



Fix). 



(A.3) 



From the common distribution. Fix), marginal distributions 
Fiixi) can be derived, as well as the mean value E[h] of quantity 
hiX): 



and compared with real data ( Krolik et al. 199^ lKaraslll997 : 



Nowak et al.|[l999l IVaughan & lJttlev..2008.) . but this would go Ftixd = f^f ]~[ dx] , E[h] = fhix) f\ f]djt:j| . (A.4) 

beyond the scope of our present work. ^, jti ^^ j=i 
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Most important for applications are the first and the second 
moments, which receive their own designation: yu,- = E[X,], and 
Rjj - E[XiXj], respectively. The meaning of Rij is the corre- 
lation matrix. Finally, it is useful to introduce the covariance 
matrix. 



where 



Qj^E[XiXj]-E[Xi]E[Xj]^Rij 



miij. 



(A.5) 



It can be proven that statistically independent random values 
are always uncorrected, in which case C,j is a diagonal matrix. 

A.2. Random processes 

Random processes are a mathematical description of measure- 
ments of a physical quantity which is either disturbed by noise 
or subject to non-deterministic evolution by itself. A random 
process is a collection of random values {X(t) : f e T) with re- 
spect to (Q, Z, P), where T is a time set (usually a subset of R or 
a set of integers). For every finite set of points {f i , . . . , f„) c T, 
one can find the distribution function /^/|,...,(„(-f i, ■ ■ ■ , x„) corre- 
sponding to the random vector X„ = iX(ti), . . ., X(t„)). 

The random process can be interpreted as a function of two 
variables, X = X{t, m). For a fixed value of some nr' e Q, 
X{t, m') is a function of time, called the trajectory or the real- 
isation of the random process. For a fixed value of f € T, the 
function X(t', m) is a real random value. A process is called sta- 
tionary if /^/,,...,r„(A:i, . . . , x„) - /^f|+r,...,f„+r(-^i, ■ ■ ■ , Xn) for every 
r e T. This implies that all moments of X(f) are independent of 
time, 

E[x\t)\ = E[xHt + r)] = Mk, (A.6) 

for all r and t. A weaker form of this condition is often used: a 
random process is called a weak-sense stationary if 

e[x(0] = E[x{t + r)], e[x\0] = E[x\f + r)]. (A.7) 

A stationary random process does not change its nature 
with time. However, in general there is no way to calculate the 
statistical characteristics of such a process knowing only a sin- 
gle realisation. For this purpose, a stronger assumption has to 
be made: a random process is called an ergodic one if for all 
fixed r e T and a real function h{x) 



1 
E\h{X{r))\^ limi {h{X{t))At. 



(A.8) 



It is commonly assumed that real processes satisfy this condi- 
tion. 

Stationary processes can be characterised by autocorrela- 
tion R{t) and autocovariance C{t) functions, 

R{t) = E[x{r) X(r + t)], C(t) = R{t) - /?, (A.9) 

which, according to stationarity, are independent of r. By set- 
ting r = -f we find that R(t) - R(-t). 

Another way of characterising processes is by their spectral 
properties, which are also our main interest in this paper Power 
spectral function of a stationary stochastic process X(t) is 



Sicj)= lim-^E[rr[^«]Ml'l, 



Tt 



1 

[x(t)] = J 



X(t) e-'"' At 



(A.ll) 



is the incomplete Fourier transform. According to the Wiener- 
Khinchin theorem the power spectral function can be calcu- 
lated from the autocovariance: Sico) - T [C(0] (w). 

A.3. Random point processes 

The concept of point processes was originally de veloped to de- 
scribe random configurations of points in space dCox & Millei 
1965l:lDaley & Vere-Jonesl2003h . One way to characterise such 
random configurations in some topological space A" c R" is by 
means of the counting measure, N(A). For every A c X, the 
counting measure gives the number of points lying in A. 

Similar to random processes, a random point process can be 
characterised by its mean value and moments. The first-order 
moment is called the intensity measure, 

Mi(A) = e[a?(A)]. (A.12) 

For every A c X, Mi (A) is the mean number of points in A. 
The second-order moment measure is then 

M2(AxB)^e[N(A)N(B)]. (A.13) 

Let {xi} be one possible configuration of points, i.e. the 
support of some A^(.). For the functions h(x) and f;(x, y) 
on X and X^, respecti vely, it follows (ICampbelll Il909 : 
Daley & Vere-Jonesll2003b 



(A. 14) 



(A. 15) 



e[^/z(;c,) = fh(x)Mi(dx), 

E\Y^g{xi,xi) = I g(x,y)M2(dxxdy). 

The concept of point process can be further generalised 
by adding marks, ki, from the mark set TC to each coordinate 
Xj from {xj]. Marks carry additional information (for example, 
they can be employed to describe the radial distribution of spots 
and some relations among spots located at diff'erent radii; see 
Sect. l3.2l i. The resulting point process on the set A'x'TCis called 
the marked point process if for every A c A" it fulfils the condi- 
tion N.{A) = N{A X TC) < oo. 

The random measure Ng{A) represents the ground process 
of the marked process A^. If the dynamics of the process is gov- 
erned only by the ground process and the marks are mutually 
independent and random values with the distribution functions 
GidK), then the measure of the marked process fulfils 



MiidxxdK) = Mgi(dx) G(d/f), 



(A.16) 



M2(dxiXd/ciXd.^;2Xd/C2) = Mg2(djciXdx2)G(dA-i)G(dA-2XA.17) 

Finally, it is useful to introduce the factorial measures, 
which will be used later to simplify various expressions. They 
satisfy definitional relations 



(A. 10) 



M[i](A) = Mi(A), 
M[2](A X B) = M2(A X B) ■ 



Mi(AnB). 



(A.18) 
(A.19) 
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